/* VanDerMonde.h
 *
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 *
 */

#ifndef _VanDerMonde_test_
#define _VanDerMonde_test_
#include <src/miscellaneous/OS_basic.h>

/*! \brief Van Der Monde solver

	It is a templated version of the numerical recipe Van Der Monde solver.
	Van der Monde system is

	\f[ \beta y = x \f]
	where \f$ x \f$ is a vector, and $\beta$ is the matrix
	\f[ \beta = \{ (\beta_i)^j \}_{i,j}  \f]
	where \f$ \{ \beta_i\}_{i} \f$ is any vector.

	To do: rewrite it in a templated iterator version for performance issues.
*/
template <class data>
data VanDerMonde(const data & x, const data & beta)
{
	OS_size i,j,k;
	data::value_type b,s,t,xx;
	OS_size size_ = beta.size();
	OS_DYNAMIC_CHECK(size_ == x.size(), "bad input in Van Der Monde system (VanDerMonde)");

	data c(size_);
	c = 0.;
	data result(c);

		c[size_-1]=-x[0];
		for (i=1;i<size_;i++)
		{
			xx=-x[i];
			for (j=size_-1-i;j<size_-1;j++)
			{
				c[j]+=xx*c[j+1];
			}
			c[size_-1]+=xx;
		}
		for (i=0;i<size_;i++)
		{
			xx=x[i];
			t=b=1;
			s=beta[size_-1];
			for (k=size_-1;k>=1;k--)
			{
				b=c[k]+xx*b;
				s+=beta[k-1]*b;
				t=xx*t+b;
			};
			OS_DYNAMIC_CHECK( t!= 0, "Singular input in Van Der Monde system (VanDerMonde)");
			result[i]=s/t;
		}
	return result;
};

#endif
